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We consider the unidirectional particle transport in a suspension of colloidal particles which 
interact with each other via a pair potential having a hard-core repulsion plus an attractive tail. 
The colloids are confined within a long narrow channel and are driven along by a DC or an AC 
external potential. In addition, the walls of the channel interact with the particles via a ratchet-like 
periodic potential. We use dynamical density functional theory to compute the average particle 
current. In the case of DC drive, we show that as the attraction strength between the colloids is 
increased beyond a critical value, the stationary density distribution of the particles loses its stability 
leading to depinning and a time dependent density profile. Attraction induced symmetry breaking 
gives rise to the coexistence of stable stationary density profiles with different spatial periods and 
time-periodic density profiles, each characterized by different values for the particle current. 

PACS numbers: 05.40.-a, 05.60.-k, 68.43.Hn 
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I. INTRODUCTION 

Much attention has been given to studying the trans- 
port of particles along narrow channels p|. Such strong 
confinement occurs for example in the ion channels of 
biological membranes Q, in zeolites and other porous 
materials @ and in microfluidic devices Q. Experimen- 
tal studies of colloidal particles confined within grooves 
etched on a surface Q have already addressed the case 
when confinement is so extreme that particles cannot 
pass one another and so single-file diffusion sets in Q. 
In many cases, the motion of a file occurs in the pres- 
ence of a periodic pinning potential. The latter can be 
induced, for instance, by defects, such as in the case of su- 
perconducting vortices moving in easy- flow channels Q, 
by other particles, as in the case of a fluctuating quasi- 1- 
dimensional (ID) channel Q, by a periodic distribution 
of charges, as in the case of the motor proteins moving 
along a microtubulus [§], or, more generally, by a peri- 
odic corrugation of the channel walls [T3, El • When the 
left-right symmetry of the pinning potential is broken, an 
externally applied center-symmetric AC drive induces a 
net drift of the file in a certain direction. The efficiency of 
such a rectification mechanism strongly depends on the 
number of particles in the file, their size, and frequency 
of the drive [l2| . 

In addition to the interaction with the channel, the col- 
loids have an excluded volume interaction between them 
[H, HH and may also exhibit mutual attraction, either 
due to Van der Waals forces [l4| or because of the pres- 
ence of other passive molecules in the solution, such as in 
the case of colloid-polymer mixtures [l5l - fl7| . It has been 
recognized that attractive forces lead to the formation of 
particle clusters and, consequentl y, to a dramatic increase 
in particle diffusion and mobility [18l ] . Such enhancement 
is explained by the mismatch between the size of the par- 
ticle clusters and the characteristic length scale of the 



corrugated potential induced by the walls of the channel. 
In the case of the diffusion of long alkane chains in zeo- 
lites, a similar phenomenon is called the "window effect" 
namely, the mobility of the alkane chain becomes 
enhanced, whenever its length is not commensurate with 
the zeolite cage. More generally, the incommensurability 
between the lateral dimensions of a biological molecule 
and the size of a catalyst is known as "shape selectivity" 
[2fj| , a recurrent scheme utilized by nature to control en- 
zymatic reactions in living cells. 

Recently, we have developed a theory [2l| , based on dy- 
namical density functional theory (DDFT) [H|-|25j], that 
captures the essential features of the condensation pro- 
cess from the disordered to the condensed state in a ran- 
domly distributed single-file of interacting particles. Pair 
attraction can be used to enhance the transport of col- 
loidal particles in ID. For instance, entrained attracting 
particles can be effectively shuttled along an asymmetric 
corrugated channel by means of a low frequency AC field. 
Collective shuttling of entrained particles directly applies 
to the problem of diffusion of long molecular chains in 
zeolites and shape selective catalytic reactions in living 
cells. In particular, we stress that collective shuttles can 
be much more efficient than some other shuttle mech- 
anisms, as they allow one to control the rate of trans- 
port by adding (removing) a single molecule to (from) 
the molecular chain. 

The present paper is organized as follows: Our main 
goal is to analyze the effects of the pair attraction be- 
tween the particles on the rectification current of single- 
files of colloidal particles which are confined within a nar- 
row channel with corrugated walls and subjected to DC 
or AC drives. The general theoretical framework for our 
analysis of our model system, which is based on DDFT, 
is presented in Scc.[TIJ The one body density distribution 
p(x, t) of the diffusing particles, which is a function of 
position x and time t, obeys a nonlinear Fokker-Planck 
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equation and DDFT 22-25] provides a closure approxi- 
mation that allows us to solve for the dynamics of p(x, t). 
The DDFT dynamical equation for the system takes the 
form of a conserved gradient dynamics, which requires 
as input a suitable approximation for the Hclmholtz free 
energy functional for the system [22-25]. In the pres- 
ence of a DC drive, the free energy contains a potential 
energy term proportional to x that acts as a continu- 
ous external energy source, i.e., the system remains per- 
manently out of thermodynamic equilibrium. In conse- 
quence, the free energy of a single-file of interacting par- 
ticles is not necessarily a monotonically decreasing func- 
tion of time. Therefore, the system can exhibit stable 
time-periodic density profiles and currents. In particu- 
lar, in Secs lIIII and IIVI we discuss the relation of the onset 
of time-periodic density variations to the condensation 
of particles into compact clusters and the depinning of 
these clusters from the channel corrugations. 

In Sec. lIV Al we examine the dynamics of point-like par- 
ticles and we demonstrate that spontaneous symmetry 
breaking induced by attraction leads to the coexistence of 
stable time-periodic and stationary densities. Multista- 
bility of the long-time density distributions indicates that 
for the same combination of parameters in our model, the 
channel can operate in two different regimes, transport- 
ing the particles with either high or low efficiency. For 
finite-size particles the range of values of the system pa- 
rameters that allows for time-periodic density profiles, is 
much broader than for point-like particles, as explained 
in Sec. UVBl 

In Sec.|V]we discuss the low frequency rectification cur- 
rent for particles driven by an AC (square wave) drive 
through a channel with a spatially asymmetric poten- 
tial. In Sec. IV Al we show that the effect of the spatial 
asymmetry is that the rectification current may be max- 
imized, using the strength of the attraction between the 
particles as the control parameter. In Sec. IV Bl we derive 
an effective equation of motion for a condensate in the 
limit of infinitely strong attraction between the particles. 
We show that the low frequency transport efficiency can 
be increased by several orders of magnitude for strongly 
attracting particles as compared to non-interacting par- 
ticles. Finally, in Sec. I VII we close with a few concluding 
remarks. 



II. DDFT FOR INTERACTING HARD RODS IN 
A PERIODIC EXTERNAL POTENTIAL 

When finite sized colloids are confined within a long 
narrow channel, the particle motion becomes one dimen- 
sional, and the particles can be modeled as ID hard rods 
of length h. We model the dynamics of N hard rods 
using overdamped stochastic equations of motion. The 
particles move in a channel of total length S and interact 
with the channel walls via a periodic corrugated poten- 
tial U(x) with spatial period L. In a system with periodic 
boundary conditions (e.g., a circular geometry), S is an 



integer multiple of L, i.e. S = ML. The integer M 
determines the average number of particles per unit cell 
of length L, to be N/M . To ensure that the combined 
length of the N rods is smaller than the total system size, 
we require Nh < S. 

The total instantaneous potential energy of the N rods, 
which move in the periodic external potential U(x) under 
the action of a time-dependent external drive A(t), is: 

where w(x) denotes the interaction potential between a 
pair of particles i and j, where i,j = 1,...,N, which 
are located at positions xi and Xj, respectively, and are 
separated by the distance x — I X £ X 7 |. In general, the 
potential w(x) can be decomposed into two terms. One 
accounts for the attraction (denoted by subscript "at") 
between the particles and the other for the hardcore re- 
pulsion (subscript "he"), i.e., w(x) = w a t(x) + Wh c (x). 
The hard-core repulsive potential ensures that the rods 
are impenetrable, that is, 



w hc (x) 



oo, x < h 
0, x > h. 



(2) 



We assume that w a t{x) becomes negligible at distances 
x much larger than a certain effective interaction range 

'int ■ 

The overdamped dynamics of the particles is described 
by N coupled Langevin equations 



1 dxi 
T~dt 



dxi 



(3) 



where T is the temperature of the system and &(t) are 
independent Gaussian white noises with the correlation 
functions = SijS(t — t'). Henceforth we set 

the solvent friction constant T = 1. 

In the case of a circular geometry, Eqs. J3} are sup- 
plemented with periodic boundary conditions (BC) with 
the period equal to the system size S. This requires that 
all functions in Eqs. ([3]) are periodic with period S. This 
is not the case if the contribution to the interaction po- 
tential w(x) from the potential tu B t(x) is long ranged, 
as particles would interact then with themselves. How- 
ever, w(x) can be made compatible with the desired pe- 
riodic BC by taking the interaction range k n t to be suffi- 
ciently small compared with the system size. Therefore, 
throughout we impose the condition w a t(S) « 0. 

The Fokker-Planck (Smoluchowski) equation for the 
time evolution of the one body density distribution p(x, t) 
is dmi]: 



dp(x,t) 
dt 



T 



d 2 p(x,t) d 



dx 2 



dx 



p(x,t) 



dU cS (x,t) 
dx 



d_ 

dx 



dx'pW(x,x',t)jj-w(\x-x'\)^) 
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where p^(x : x',t) is the non-equilibrium two-body dis- 
tribution function for the particles in the system and 
Ueg(x,t) — U(x) — A(t)x is the effective external poten- 
tial. Note that the density profile p(x, t) is normalized, so 
that the spatial integral over p(x, t) is equal to the total 

number of particles in the system; i.e. /_5^y 2 P( x > ^) dx ~ 
N. 

In order to solve Eq. (TJJ , a suitable closure approxima- 
tion for p( 2 ) (x, x' , t) is required. The approach taken in 
DDFT is to approximate p^ (x, x' ,t) by the two-body 
distribution function of an equilibrium fluid with the 
same one-body density profile as the non-equilibrium sys- 
tem [22l - [25| . This closure relates the integral in the final 
term in Eq. Q to the functional derivative of the excess 
part of the Helmholtz free energy functional F[p], which 
is the central quantity of interest in equilibrium density 
functional theory [lj, [2(| [53] ■ Making this approxima- 
tion yields the following equation for the dynamics of the 
one-particle density distribution p(x,t): 



dp(x,t) d 
dt dx 



p(x,t) 



d SF[p(x,t)} 
dx 5p{x,i) 



(5) 



For the case with periodic BC, the Helmholtz free energy 
functional F[p] is of the form QUI]: 

F[p(x,t)] = T / dxp(x,t)[]np(x,t)-i\ 



(0) 



/ dxU e g(x,t)p(x,i) 
F hc [p] + F at [p], 



where the first term on the right hand side is the ideal-gas 
contribution to the free energy, Fh c is the excess contri- 
bution to the free energy due to the hardcore repulsion 
between the particles, and F at represents the contribu- 
tion due to the attractions between the particles. In a 
mean-field approximation [26| , F a t is given by 



fAp\ 



dx 



x+% 



dx'w at (| x — x' \)p(x)p(x'). 



(7) 

The exact expression for the equilibrium excess 
Helmholtz free energy for hard rods of length h, Fh c [p], 
was first presented in Ref. (28|. The result is: 

Fhc\p] = \J \ dx {p (* + 5) + P ~ |) ] .(8) 

where <fi[p(x,t)] = — Tin [1 — r)(x, t)] and r)(x, t) = 
Jx-h/2 p( x ' >t)- It should be emphasized that the 
functional in Eq. (|HJ|, is strictly only exact for ID equilib- 
rium systems of hard-rods treated in the grand canonical 
ensemble [26|, [28|. However, as the (average) number 
of particles in a system is increased, the difference be- 
tween results from treating a system canonically or grand 



canonically diminishes, so that the theory can safely be 
extended to describe single-files with a fixed but large 
numbers of rods. 

In earlier work, the DDFT approach was used to study 
the dynamics of an ensemble of pur e hard rods (i.e. with 
no attractive interactions) [12, l23l f29j ] . More recently 
plj . we have applied the DDFT formalism to describe 
a file of hard rods interacting via a pair potential with 
an attractive contribution. In both cases, the free en- 
ergy functional given in Eq. (|5]) was shown to reproduce 
fairly closely the results from Brownian dynamics com- 
puter simulations [i.e. results from numerically integrat- 
ing Eqs. (|3[)], even for relatively small numbers of parti- 
cles, N > 10. 

Using Eqs. © and ©, Eq. ([5]) can be rewritten in the 
form of a conservation law, i.e. in terms of the instanta- 
neous current density J(x,t), 



dp(x,t) 



dt 



dJ{x,t) 
dx 



(9) 



where 



J(x,t) = p(x,t) 



T 



1 



-T— \np(x,t) 

p{x + h, t) 
-r](x + h/2,t) ~ 



dU{x) 
dx 

p(x - 



+ A(t) 
h,t) 



dx' p(x',t) at (x 
dx 



1 — rj(x 
x>) 



h/2,t) 



(10) 



Before we discuss our results for solutions of Eqs. © 
and (fTU)) . we would like to point out that we see many 
parallels between the dynamics described by the partial 
integro-differential Eqs. (JS|) with ^ and the dynamics of 
partially wetting drops and films on solid substrates de- 
scribed by so-called thin film equations (fourth order par- 
tial differential equations) [30(. Formal similarities be- 
tween thin film equations and some Fokker-Planck equa- 
tions for interacting particles have recently been pointed 
out in the case of spatially-asymmetric ratchets with a 
temporal AC drive [31|. In the present case, we find that 
some of our results are similar to results for the case 
of liquid drops on horizontal [32| and inclined [Hj] het- 
erogenous solid substrates. These similarities arise due 
to (i) the similar gradient dynamics form of the evolu- 
tion equation for the conserved field (here p), (ii) the 
presence of similar physical effects. The role of the at- 
tractive and repulsive force between particles is taken by 
the partial wettability of the liquid [HI; the stabilizing 
role of diffusion is played by surface tension; the channel 
corrugations are similar to substrate heterogeneities; the 
DC drive corresponds to constant driving parallel to the 
substrate (e.g., drop on an incline); and the present AC 
drive is similar to substrate vibrations [35[ or oscillating 
electric fields [36j . 

Here, we solve Eq. ([9]) imposing periodic BC on the 
domain x £ [-S/2, S/2] which is centered at the origin, 
x = 0. Due to the external driving force the system re- 
mains permanently out of thermodynamic equilibrium. 
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Therefore, in an infinite domain or in a finite domain 
with periodic BC, the non-equilibrium dynamics of the 
system with DC drive is not relaxational, in contrast to 
the case of zero-flux BC that would correspond to a closed 
finite system [22J. Note, however, that the channel cor- 
rugations may result in local equilibria. An AC drive 
keeps the system out of equilibrium for any BC. The non- 
relaxational character is reflected by the finding that with 
periodic BC, the free energy F[p(x,t)] in Eq. ([6]) is not 
necessarily a monotonically decreasing function of time, 
i.e., it does not play the role of a Lyapunov functional for 
the gradient dynamics Eq. (JSJ) - This can best be seen by 
considering the total time derivative of the free energy 



dF[p] 
dt 



5F dp(x, t) 



s 5p 

2 ' 

% 5F d 

_s 5p dx 

2 ' 

5F d 5F 
Sp ^ dx Sp 



dt 



dx 



p(x,t) 



d_6£ 

dx Sp 



P 



dx 

d_SF_ 

dx bp 



efat.ll) 



where we have used Eq. ([5]) and integrated by parts. The 
functional derivative SF/5p is not periodic in x, due to 
the "tilted" effective potential U c s = U(x) — A(t)x, and 
so the boundary term in the last line of Eq. (ITT1) does not 
normally vanish. It is equal to ASJg/ 2 > 0, where J5/2 
is the current density on the boundary, and Js/2 < for 
A < 0. It then follows that for A^O the time derivative 
dF[p(x 1 t)]/dt is not necessarily a negative quantity. This 
allows for time-oscillatory behavior of the system, even in 
the long-time limit. In other words, Eq. ©, with A(t) = 
A, may admit stable cyclo-stationary or time-periodic 
solutions, the existence of which will be discussed in the 
next section. 

For stable time-periodic solutions, the average particle 
current J is obtained from Eq. ([TO]) by averaging J(x, t) 
over position and over time, namely, 



J 



1 



dx 



t'+T 



dt J(x, t), 



(12) 



where r is the period of the oscillations. Throughout 
this paper we will use the average current per particle 
J, which is related to the total current in Eq. (IT2l via 
J = J/N. 



III. SPONTANEOUS CONDENSATION 

In the absence of the channel potential and without 
AC or DC drive, i.e for U(x) cS = 0, the DDFT equation 
([5]) admits a stationary homogeneous solution p(x) = pa, 
with constant average density p — N/S. Because of the 
competition between the destabilizing attractive forces 
and the stabilizing effect of the thermal motion of the 
particles (diffusion), the homogeneous density distribu- 
tion is not always stable. In order to minimize their free 



energy, the attracting particles tend to condense into a 
set of clusters. This tendency is opposed by the stabiliz- 
ing action of diffusion, which leads to a spreading of the 
particles away from one another. Depending on which 
process dominates, the homogeneous state p(x) = po may 
either be (linearly) stable or unstable. Note that in this 
case the dynamics is relaxational and the functional ([5]) 
is a Lyapunov functional [cf. Eq. Qlip] . 

One may alternatively take a thermodynamic rather 
than a dynamical point of view for understanding this 
instability: Recall that the Helmholtz free energy of the 
system F = U — TS, where U = (<&) is the internal en- 
ergy, S is the entropy of the system and (• • • ) denotes a 
statistical average [14j | . At equilibrium, F is minimal. At 
high temperatures T, this is achieved by maximizing S 
(i.e. by dispersing the particles throughout the system in 
a maximally disordered way), since the term — TS is the 
dominant contribution to the free energy at high tem- 
peratures T. However, at low temperatures, the internal 
energy contribution hi dominates the free energy F and 
so the attracting particles minimize the free energy by 
gathering together to minimize the internal energy IA. 

We model the attractive contribution to the pair poten- 
tial between the particles by a simple exponential func- 
tion of the form 



Wa. t (x) = -aexp(-Ax), 



(13) 



where the parameters a and A characterize the strength 
of attraction and the attraction range Z a t = 1/A, respec- 
tively. It should be noted that the mean field approx- 
imation for the contribution to the Helmholtz free en- 
ergy due to the attractive interactions (J7J) used here, is 
only quantitatively reliable when Xh < 1, i.e. when any 
given particle is interacting with several of its neighbors 
so that a mean-field approximation is appropriate. How- 
ever, even when the attraction range is somewhat shorter 
than this, the mean-field approximation remains qualita- 
tively correct. 

In order to study the stability of the homogeneous 
state, we linearize Eq. ([5]) using the standard plane wave 
anstatz p(x,t) = po + ee^i k ) t +' tkx , where the sign of the 
growth rate /3(k) determines the stability of the homoge- 
neous solution po. Note that for the stable fluid the dis- 
persion relation /3(k) is closely related to the static struc- 
ture factor S(k), via f3(k) = -Tk 2 /S(k) [H[. Neglecting 
exponentially small terms of the order of exp (—AS), we 
obtain to leading order in e 



-Tk l 



2T 'kp sin (kh) 4Tpg[sin (kh/2)] 



1 - Poh 



(I-Pohy 



2a\pok 2 
k 2 + A 2 



(14) 



Inspection of Eq. (TT4")) shows that two different scenar- 
ios exist where the homogeneous density po is linearly 
unstable. On the one hand, the solution p(x, t) = p can 
become unstable via the standard spinodal phase separa- 
tion mechanism where the material separates into regions 
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FIG. 1: (Color online) (a) Typical dispersion relations /3(fc) 
close to the onset of the spinodal mode (sp), as shown by 
dashed lines and the freezing mode (fr), as shown by solid 
lines. The /3(fc) curves near the spinodal instability are for 
po = 0.2 and a = 10, 11.8 and 17, which are close to the 
onset of the instability. The curves for the freezing instability 
are for po = 0.6 and a — 4, 4.9 and 6. The other parameters 
are T = 1, A = 3, h = 1. (b-c) Stability diagram for the 
system with uniform constant density p(x) = po, in the plane 
spanned by T = ThA/a and p = poh for (b) £ = Ah = 3 
and (c) £ = 1.5. The labels "sp" and "fr" denote (hatched) 
regions where the uniform density is only linearly unstable 
to the spinodal and freezing instability, respectively. In the 
cross-hatched region both modes are unstable. 



of low density (gas) and hight density (liquid) , which we 
call the "spinodal mode" (denoted by "sp"). It is associ- 
ated with an instability against harmonic perturbations 
with wave numbers < k < fco, where kg is an upper 
value obtained by solving the equation /3(fco) = 0. The 
spinodal instability sets in when the leading coefficient 
of the expansion of f3(k) in powers of k 2 (i.e., the coeffi- 
cient of the k 2 term) vanishes, i.e., the wave number at 
onset is zero. This corresponds to a type II S instability 
in the classification of Cross and Hohenberg [37j ■ Based 
on the Taylor expansion of the relation (jT4")l one obtains 
T = 2/5(1 — p 2 ) as the condition for the onset of the 
spinodal mode [cf. Fig. H{b)]. Here, we have introduced 
the reduced temperature T — ThX/a and the reduced 
density p — p^h. The critical point for decomposition is 
found at f c = 8/27 and p c = 1/3. 

The stability condition is equivalent to the thermo- 
dynamic stability criterion requiring that the isother- 
mal compressibility be negative. This corresponds to 
the Helmholtz free energy per unit length of the system 
becoming concave, namely, the boundary of the spin- 
odal instability is given by S 2 F[p]/5p 2 \ Po — 0. When 
5 2 F[p]/Sp 2 \ Po > the uniform solution is linearly sta- 
ble and when S 2 F[p]/Sp 2 \ Po < the uniform solution is 



linearly unstable. The dispersion relation f3(k) close to 
the onset of the spinodal mode is shown in Figfjja) by 
dashed lines, which correspond to three different values 
of the attraction strength a, chosen close to the spinodal 
instability threshold as obtained from Eq. (IT4l) . 



On the other hand, the solution p{x 1 1) = po can be- 
come unstable via a freezing mode of the system, where 
the particles become localized and the density profile 
exhibits a series of sharp peaks separated by distances 
smaller than L. This instability corresponds to type / s in 
the classification of Ref . [37| • In the context of reaction- 
diffusion systems it is sometimes referred to as a Turing 
instability 38,]. This mode (which we denote by "fr") 
sets in at a non-zero critical wave number fc c , as illus- 
trated in FigOJa) by the solid lines for po = 0.6. Beyond 
onset, the freezing mode gives rise to the growth of peri- 
odic modulations in the density profile with a wavelength 
w 2n/k c . 

In Figs.QJb) and (c) we plot the linear stability dia- 
gram of the system with uniform density po, in the plane 
spanned by T — ThX/a and p = poh for the interaction 
length ratio (b) £ = Xh = 3 and (c) £ — 1.5. The labels 
"sp" and "fr" mark regions where the uniform system is 
linearly unstable to the spinodal and freezing mode, re- 
spectively. As discussed above, the onset of the spinodal 
mode only depends on the reduced temperature T and 
reduced density p. In contrast, the onset of the freez- 
ing mode also depends on the value of £. For relatively 
long rods with £ = 3, where the attraction range is short 
compared to the core size h, the region of the freezing 
instability extends down to moderate values of the re- 
duced density p 0.5. For a reduced temperature above 
T c = 8/27 only the freezing instability exists at large and 
moderate p. For small values of £, corresponding to the 
attraction range A -1 being large compared to h, and for 
which the mean-field approximation for the free energy 
used is expected to be most reliable, the freezing mode 
is only found at extremely high packing fractions, p « 1. 
At such high densities, the critical wave number k c of the 
freezing mode is approximately k c w 2ir/h, giving rise to 
the formation of density peaks separated by a distance 
~ h, as one would expect for a frozen system. 

Below T c , the spinodal mode exists for a range of p that 
with decreasing T extends on both sides of the critical 
value p — 1/3. This implies that the spinodal mode sets 
in at smaller and smaller values of the density as T is de- 
creased. At large p there exists a region where both linear 
modes are unstable [cross-hatched in Figs.[TJb) and (c)]. 
The region where only the freezing mode exists is shifted 
towards higher densities as T decreases. The definition 
of T implies that for any given physical temperature T, a 
decrease in the interaction strength a below the thresh- 
old value a c = 27THX/8 stabilizes the spinodal mode. 
For any finite temperature the system can be quenched 
into the freezing unstable region by increasing the aver- 
age density in the system. 
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IV. DC DRIVE 

Having considered the stability of the uniform system, 
we now consider the non-uniform system that is subject 
to a periodic external potential U (x) and the DC driving 
force A. We model the periodic potential, induced by the 
corrugated channel walls, by the standard bi-harmonic 
ratchet potential U(x) = sin (2ttx) + 0.25 sin (4ttx) 
Recall that to drag a single particle over one of the barri- 
ers in U(x), one must apply a force Ar = 3ir to pull 
the particle over the barrier to the right and a force 
Al = 37r/2 to pull it to the left. Al and Ar are 
termed the left and the right depinning thresholds, re- 
spectively Note that all the results reported in this 
section for a DC drive remain qualitatively valid even for 
a simpler symmetric periodic external potential, such as 
U(x) = sin (2irx). The case of a periodic external poten- 
tial without and with DC drive shows some similarities 
to liquid drops/films on periodically heterogeneous sub- 
strates without [32| and with [39[ a driving force parallel 
to the substrate, respectively. The former case will be 
explored elsewhere. Below we discuss similarities and 
differences for the case with DC driving. 



A. Zero rod length and finite interaction range 

As a reference system we first consider a file of point- 
like particles, i.e., with h = 0, interacting solely via the 
exponential soft core potential H7 a t(x) and driven by a 
DC external force. When the characteristic interaction 
range l a t between the point-like rods is small, i.e. when 
lax = 1/A — > 0, a local approximation can be made for 
the dynamical equations for the system, as shown in 
Refs. [4(| |4l|. In this limit, the integral involving to at 
in Eq. (fTU)) can be reduced to a local function of the form 
~ gp{x,t)dp{x,t) I 'dx, where the coefficient g is a pa- 
rameter determined by the str eng th of the interactions 
between the particles. In Refs. [40|, |4l| it was shown that 
when g is increased beyond a certain critical value, the 
density distribution of the attracting particles exhibits 
a spontaneous symmetry breaking transition, where the 
stationary periodic density profile with period L [the pe- 
riod of the modulations in U(x)] becomes unstable and 
evolves toward a stable stationary distribution with pe- 
riod S (the total sy_stem length) . We now go beyond the 
analysis of Refs. [40t l4lj and consider such a symmetry 
breaking mechanism for particles interacting via a poten- 
tial with a nonzero interaction range, Z at 7^ 0. 

For convenience we fix the average particle density to 
be pL = 1, corresponding to one particle per period L of 
the channel, and we set the constant drive, A = — 1. Us- 
ing the numerical continuation package AUTO we 
follow the branch of solutions corresponding to a station- 
ary density distribution p s (x), that originates from the 
stationary density profile for the case when a = (i.e. a 
non- interacting ideal-gas of particles). Note that in the 
long-time limit, the density profile for the ideal-gas re- 



mains stable and stationary, regardless of the form of the 
channel potential, U(x), the magnitude of the drive, A, 
or the temperature of the system, T [H, Q • 

We determine the stationary solutions of Eq. (J9)) with 
the current given by Eq. QlOjl that contains nonlocal 
terms, by means of the Fourier mode method described 
in Ref. [45J. The density profile is discretized over the 
domain [—5/2, 5/2], derivatives are obtained using fi- 
nite difference approximations, and the non-local terms 
are calculated using a Fast Fourier transform. We start 
from the equation for the stationary solution of Eq. © , 
dJ(x,t)/dx = 0, which is then written as a set of alge- 
braic equations for the Fourier components of the cur- 
rent J(x, t) and from this we obtain our solutions for the 
stationary density profile p s {x). Using the continuation 
package A UTO allows us to detect the presence of Hopf 
bifurcations as well as to trace the solution branches for 
both the stationary solutions and the time-periodic ones 
that emerge from them. 

We begin by discussing the bifurcation diagrams of the 
stationary solutions of Eq. Q on varying the interaction 
strength, a, for a fixed value of the range parameter of 
the pair potential, A = 5, and the fixed U(x). These are 
shown in Figs.[5Ja)-(c), for three different systems with 
lengths, S = 2L, 3L and 4L, respectively. The solid lines 
correspond to stable solutions and the dashed lines to 
unstable (saddle point) solutions. The labels "HB" and 
"BP" stand for Hopf bifurcation and branching point, 
respectively. 

As the interaction strength is increased beyond a crit- 
ical value ct c , the lL-periodic solution that is stable for 
small a becomes unstable either via a (period-doubling) 
pitchfork bifurcation (for S = 2L), or via a Hopf bifurca- 
tion (for S > 2L). In the case of the pitchfork bifurcation, 
displayed in Fig. 0(a), a double branch of stable solutions 
emerges at the bifurcation point. The two branches are 
related by the discrete translation symmetry x — > x + L 
and can therefore not be distinguished in Fig.[3Ja). This 
new branch corresponds to a solution with a larger spatial 
period, equal to the system size S = 2L, and a smaller 
value for the particle current J. One may say that for 
a > a c , the periodic potential is not strong enough to pin 
the clusters against their natural tendency to coarsen. 
For S — 2L, we know on general grounds that [in a ho- 
mogeneous system without driving A = U(x) = 0] there 
are two possible coarsening modes: a translation mode, 
where the two clusters move towards each other, and a 
volume transfer mode, where material is transfered from 
one cluster to the other 30]. It is known that both are 
stabilized by substrate heterogeneities exerting a strong 
enough pinning influence (32| . Not much is known, how- 
ever, for driven systems (A ^ 0). Our DDFT simulations 
show that in the present system, the instability is related 
to the volume transfer mode of coarsening. The corre- 
sponding stable (solid line) and unstable (dashed line) 
stationary density profiles are displayed in Fig.[2jd) for 
the case when a — 2.2. 

The bifurcation diagram for the system with length 
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FIG. 2: (Color online) Panels (a), (b) and (c) display bifurcation diagrams for the stationary density distributions in terms of 
the average current J versus the interaction strength a, for A = 5, p = 1, A = —1, h = and total system length (a) S = 2L, (b) 
3L and (c) 4L. The solid and dashed lines correspond to stable and unstable solutions, respectively. The points labeled "HB" 
and "BP" denote Hopf bifurcation and branching points, respectively. Panels (d), (e) and (f) display selected corresponding 
density profiles. The solid lines are profiles with spatial periods 2L, 3L and 4L, respectively, for a = 2.2. The dashed lines 
represent the unstable solution with period L. In (f) the solid and dotted lines represent all four possible 4L-periodic solutions. 



S = 3L is qualitatively different from the one for the 
S = 2L case, as can be seen in Fig.Efb). One ob- 
serves that the lL-periodic solution becomes unstable via 
a Hopf bifurcation. There exist branches of stationary so- 
lutions where the x — > x + L translational symmetry is 
broken. However, they do not touch the primary branch 
of the li-periodic solutions but are generated through 
a saddle-node bifurcation at the point marked by "LP" . 
Solutions on these branches have period S = 3L and are 
either stable (upper branch) or unstable (lower branch). 
An example of a stable 3L-periodic solution for a = 2.2 
is displayed in Fig.^e). 

The bifurcation scenario for S = AL, displayed in 
Fig.[2jc), is substantially more complex. The lX-periodic 
solution becomes unstable through a Hopf bifurcation. 
Very close to the HB point, the unstable stationary so- 
lution undergoes a primary period-doubling pitchfork bi- 
furcation BP. Note that this first BP point for S = AL 
coincides as expected with the BP point for S = 2L. The 
newly formed 2L-periodic solution is unstable and under- 
goes a further period-doubling pitchfork bifurcation at a 
second BP, which lies very close to the first BP, as shown 
in the inset of Fig.^c) . 

For interaction strengths significantly larger than the 
critical values corresponding to the BP and the HB 
points, the only stable stationary solution of the DDFT 
equation ((SJ) has a period equal to the total system size S. 



Note that the multiplicity of the branch of the solutions 
with broken x — > x + L symmetry depends on the total 
system size. For instance, for S = AL there exist four 
such branches with spatial period S 1 , as can be seen in 
Fig.[2jf). Each solution exhibits a prominent maximum 
centered around one of the four minima of the external 
potential U(x). The solutions on the 4 branches are re- 
lated by the symmetry x — > x + L. Therefore all of them 
correspond to the same value for the particle current J 
and they can not be distinguished from one another in 
Fig.^c). 

From the results displayed in Fig. [5] we may draw two 
important conclusions: First, the detected Hopf bifurca- 
tions of the pinned stationary solutions signals the onset 
of time-periodic solutions of the DDFT equation ([S]) , even 
in the presence of a time independent drive. Second, for 
certain values of the interaction strength a, two stable 
stationary solutions may coexist, giving rise to current 
multiplicity. 

These two findings are illustrated in detail in Fig.|3l 
where we display a magnification of the region close to 
the bifurcations in Fig.[^b). In addition to the Hopf 
bifurcation (HB) and the saddle-node bifurcation (LP) 
of the stationary 3L-periodic solutions, we display the 
branch of time-periodic solutions of Eq. ([5]). It emerges 
at the HB point and terminates in a homoclinic bifurca- 
tion (labeled by "hm") where the time-periodic solution 
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FIG. 3: (Color online) Magnification of the region of Fig.[2jb) 
close to HB point. The line connecting HB and "hm" corre- 
sponds to stable time-periodic solutions of Eq. ©. The label 
"hm" stands for homoclinic bifurcation point. (Inset) Tem- 
poral period r of the time-periodic solutions as a function of 
aiim — a. 



(limit cycle) collides simultaneously with all three unsta- 
ble 3-L-periodic solution (unstable equilibria) 46]. The 
inset of Fig. [3] gives the temporal period r as a function 
of the distance to the homoclinic bifurcation a^m — a. It 
shows a logarithmic dependence as expected close to a 
homoclinic bifurcation. We emphasize that these time- 
periodic solutions are stable, i.e., the corresponding Flo- 
quet multipliers are always located within the unit circle 
(not shown). Note that for clarity we not only suppress 
the branch of time-periodic solutions in Fig.[5Jb) but also 
a similar branch in Fig.^c), for the system with S = 4L. 

In nonequilibrium driven systems, the loss of stability 
of the stationary solutions and the appearance of time- 
periodic solutions with a larger mean flow is sometimes 
associated with the concept of depinning. For example, in 
the study of liquid droplets on an inclined heterogeneous 
solid substrate, the dynamics of drop depinning has been 
studied in great detail - see for example Refs. [3!j|47[ and 
references therein. In this situation the depinning is gen- 
erally a transition from a steady droplet, pinned by the 
heterogeneity of the substrate, to a moving droplet, slid- 
ing down the incline under the action of gravity (or other 
driving forces parallel to the substrate) . The depinning is 
usually investigated by increasing the driving force with 
all other parameters kept fixed. In such a case the dom- 
inant depinning mechanism is often related to a Saddle 
Node Infinite PERiod (sniper) bifurcation, although de- 
pinning via a Hopf bifurcation may also be observed in 
certain parameter regions [H, [3f| [43] . The depinning ex- 
hibited by the present system is observed when increasing 
the particle attraction a, for a fixed value of the exter- 
nal drive A and potential U(x). This would correspond 
to a decrease in wettability for a droplet depinning in 
the thin film model. Note also that this collective depin- 
ning is very distinct from the T = transition that is 




FIG. 4: (Color online) (a) Snapshots of the time-periodic den- 
sity, p(x,t), for a = 2.08 at three times, to (dashed line), 
to + r/3 (dot-dashed line), and to + 2r/3 (solid line). The tem- 
poral period of the solution is r = 30; the remaining numerical 
parameters are as in Fig.[2jb). The channel potential U(x) is 
also displayed (heavy solid line limiting the shaded area), (b) 
The time dependent current J(t) = (1/N) J^ s ^ 2 J(x,t) dx, 
corresponding to the solution in (a). 



also referred to as 'depinning', when the drive on a single 
particle exceeds either Al or Ar, the left and right single 
particle depinning thresholds 

At the HB point, the newly formed stable time- 
periodic solution has a finite period, as it can be seen 
from Fig-[3j In order to illustrate the dynamics of the 
depinning of the stationary solution, we set a = 2.08 
and plot in Fig.@|a) snapshots of the time-periodic solu- 
tion p(x,t) at three subsequent times, to, to + t/3, and 
to + 2t/3, where to was chosen as described below and r 
is the temporal period of the solution. Inspection of the 
density profiles indicates that the depinned solution can 
be seen as a superposition of two parts: a stationary part 
with spatial period L and a time-periodic part with spa- 
tial period S = 3L, that slides 'on top' of the stationary 
part. The time-periodic part corresponds to a wave trav- 
eling to the left, that is, in the direction of our negative 
constant drive. Here, at t — to, the absolute maximum of 
the density profile is located at the rightmost minimum 
of the channel potential, U(x). After one third of the 
temporal period t, the absolute maximum has moved to 
the central well of the channel and after two thirds of r, 
the maximum has finally reached the leftmost well. After 
one full period r, the cycle is repeated. 

The time-periodic solution changes its character along 
the branch in a continuous manner. With increasing at- 
traction strength the amplitude of the time-periodic part 
becomes larger as compared to the steady part until fi- 
nally most of the particles travel. They travel, however, 
not in the form of a translation of a compact cluster, 
but rather in the form of a volume transfer of the clus- 
ter from one potential well to the next. The temporal 
period becomes larger with increasing a and the overall 
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flux oscillates between a low absolute value (when the 
cluster sits in a well) and a large absolute value (when 
the cluster is transferred to the next well) . This is shown 
in Fig. 0|b). With increasing a the dependence of the 
flux on time becomes increasingly non-harmonic as the 
cluster spends an increasing fraction of the time period 
around the three maxima of the channel potential. In the 
vicinity of the homoclinic bifurcation the density profiles 
for clusters mainly localised at one of the three maxima 
closely resemble the corresponding profiles on the three 
unstable stationary 3L-periodic solutions. This also im- 
plies that at the homoclinic bifurcation the stable cycle 
collides with all three unstable equilibria at once. 

Summarizing the results displayed in Figs. and [H 
we conclude that, for values of a between the points 
labeled by LP and HB, there exist two stable station- 
ary solutions, with spatial period L and 3L, respectively. 
Moreover, between the points HB and hm, a stable time- 
periodic solution coexists with the stable stationary 3L- 
periodic solutions. By perturbing the time-periodic den- 
sity profile with a finite amplitude disturbance, one can 
induce the transition to the stable stationary 3L-periodic 
solution. To do so, one starts a simulation in time with 
a stable time-periodic solution and adds a finite (mass- 
conserving) perturbation. If the perturbation is large 
enough, the solution evolves after a short transient to- 
ward the stable 3L periodic stationary solution. Note 
also that we were not able to find the opposite transition: 
Perturbing the stable stationary 3L-periodic solution by 
shifting it slightly in the direction of the drive will 'depin' 
the cluster only for a short transient. It moves to the left 
and settles into the next potential well, i.e., it moves to 
the stable stationary 3L-periodic branch that is related 
by the translation x — > x — L. 



B. Finite rod length and short range attraction 

We discuss now the effects of having a finite rod length 
in addition to the attraction between the particles. In 
Fig.[5ja) we display the bifurcation diagram in terms of 
the stationary current J as a function of a for rod lengths 
h = 0,0.1, and 0.2 for a domain length S = AL. A rel- 
atively small change in the size of the rods is sufficient 
to cause a significant change in the current. First, one 
observes that the magnitude of the current at a = 
increases with h; this also remains true for a > 0. Sec- 
ond, the critical interaction strength, a c , at which the 
lL-periodic solution looses its stability via a Hopf bifur- 
cation, increases with h; i.e. as expected a system of finite 
length rods is more stable than the reference system with 
h — > 0. Third, the region in parameter space in which 
the stationary lL-periodic and the 4L-periodic solutions 
coexist, shrinks as h is increased. 

This can be explained as follows: For the 4L-periodic 
solution to be stable at relatively small values of a, one 
must squeeze all the particles (there are 4 particles in the 
system with length S — AL and p — 1) into a small part 



of the total system, not larger than half a ratchet period, 
Ah < L/2. As a consequence, the critical rod length 
above which the 4L-periodic and the lL-periodic solu- 
tions are unlikely to coexist, is approximately h = 0.125, 
for L = 1. A magnification of the bifurcation diagram for 
h = 0.1 slightly below this critical value is displayed in 
Fig. 03b). There, the stationary lL-periodic solutions be- 
come unstable at the Hopf bifurcation (HB) and a stable 
branch of time-periodic density profiles emerges super- 
critically [heavy green solid line in Fig.[5^b)] . Slightly 
beyond the Hopf bifurcation, the unstable branch of sta- 
tionary li-periodic solutions undergoes a supercritical 
period-doubling pitchfork bifurcation (BP) (a « 3.35). 
The emerging branch of stationary 2L-periodic solutions 
is unstable w.r.t. two modes. It becomes more unstable 
at a secondary period-doubling pitchfork bifurcation at 
(a Ri 3.45) (BP). The bifurcating branch consists of sta- 
tionary unstable 4L-periodic solutions with 2 unstable 
eigenmodes. One of them is stabilized at a first saddle- 
node bifurcation at a ~ 3.55 where the branch turns back 
toward smaller a. The branch of stationary 4L-periodic 
solutions finally becomes stable at another saddle-node 
bifurcation at a w 3.47, where it turns again towards 
larger a. The branch of stable time-periodic solutions 
terminates as in the case of h = length rods in a ho- 
moclinic bifurcation on the branch of unstable stationary 
4L-periodic solutions. The exact location of the homo- 
clinic bifurcation (labeled "hm") is very close to (but 
numerically clearly distinguished from) the saddle-node 
bifurcation. The temporal period of the time-periodic so- 
lutions diverges logorithmically on approaching the "hm" 
point, as shown in Fig.[5jf). 

To obtain some indication as to which stable solution 
might be selected in time evolutions of the DDFT, start- 
ing from various initial states, we compute the (time- 
averaged) Helmholtz free energy per unit length, (F/ S) , 
for all stable solutions. They are displayed in Fig.[5Jc) 
and (d) for h = 0.2 and h = 0.1, respectively. In calculat- 
ing these, we subtract the non-periodic potential energy 
term, f_ S / 2 Ap(x) dx, associated with the DC drive, from 
the full expression in Eq. ©. This ensures that solutions 
on branches that are related by the discrete x — > x + L 
translational symmetry have an identical value for the 
free energy under periodic BC. 

For h — 0.2, we observe in Fig.[SJa) that there exists 
no branch of stationary 4i-pcriodic solutions; instead the 
stable branch of 2L-periodic solutions continues toward 
large a. Note that this branch is unstable when it bi- 
furcates from the 1L solutions, but becomes stable as 
a result of a another Hopf bifurcation at a ss 5.6. The 
emerging time-periodic branch is unstable and will not be 
further considered here. The only stable solutions with 
spatial period equal to the system size, S = AL, are the 
time-periodic ones, which correspond along most of the 
branch to a single compact cluster of particles traveling in 
the direction of the drive. Close to the Hopf bifurcation, 
it resembles a small amplitude wave moving 'on top' of 
the stationary 1L state. Further away from the bifurca- 
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FIG. 5: (Color online) (a) Current J versus a for h = 0, 0.1 and 0.2, as indicated on the curves. The other system parameters 
are T = 1, A = 5, S = 4L, A — —1 and p = 1. The symbols denote the Hopf bifurcation (HB) and the branching points (BP). 
The primary bifurcation is always a HB. The solid and dotted lines for h — 0.1 and h — 0.2 represent stable and unstable 
solution branches, respectively. For h = 0.2, the branch originating from the HB point is the branch of stable time-periodic 
solutions. Next to each stable branch is a label indicating the spatial periodicity of the corresponding solutions, (b) displays 
a magnification of the region in the vicinity of the HB point for h — 0.1. The branch of stable time-periodic solutions starts 
at the HB point and terminates at the hm point. Panel (d) shows the average free energy (F/S) for h = 0.1 and the same 
range of values of a as in (b). Panel (c) shows the free energy of the various solution branches, for h = 0.2. Panel (e) shows 
a snapshot of a typical time-periodic solution, obtained for parameters as in (c) and a = 8. Panel (f) represents the temporal 
period of the stable time-periodic solutions in (b) as a function of the distance from the hm point, i.e. (oa lm — a). Dashed line 
is the decay law r ~ In (Qhm — a)- 



tion the behaviour resembles the one described above in 
connection with Fig. 2) Most of the particles travel in the 
form of a volume transfer of the cluster from one potential 
well to the next. Increasing a further, at about a « 7 the 
flux increases by about 50% over a very small cv-range. 
And the cluster morphology also changes from a compact 
"drop-like" shape to a multi-hump localized structure as 
depicted in Fig.[5je) , with an arrow indicating the direc- 
tion of motion of the cluster. Each hump corresponds to 
a single particle. The particles in the cluster are strongly 
bound together and the distance between the particles 
remains almost constant as the cluster moves through 
the system as a single unit. This implies that at a ~ 7 
the transport mode also changes from a volume transfer 
mode or to a translation mode. 

In contrast to the case h — 0.1, the time-periodic 
branch continues toward large a. In other words, for 
h = 0.2 all 4L-periodic solutions are depinned. In 
Fig-EJc) we see that time-periodic solutions have on av- 
erage a lower free energy than the stable stationary 2L- 
periodic solutions. However, as the system is perma- 
nently out of equilibrium, in general, the solution of lower 



free energy is not necessarily the one that the system 
converges to in the long time t — > oo limit. Thus, for 
h = 0.2 the onset of the time-periodic solutions of the 
DDFT equation is associated with a transition between 
two major transport modes: (i) At small values of the 
attraction strength a or, equivalently, for high temper- 
atures, stationary density distributions exist, with the 
particles uniformly distributed among the wells of the 
channel potential. Under the action of the stochastic 
(thermal) noise, the particles jump occasionally either 
to the right or to the left, but with a higher probabil- 
ity for jumps in the direction of the applied drive. One 
may call this the "stationary mode" . (ii) At larger a (or 
smaller temperature) , time-periodic density profiles seem 
to dominate. They either correspond to transport from 
well to well by a volume transfer mode or by a translation 
mode. The latter correspond to depinned compact clus- 
ters in which strongly attracting particles travel together. 
One may call this the "condensed traveling mode" . Such 
a traveling cluster has a characteristic length ~ hN and, 
in the limit where the attraction a is strong (i.e. when 
a 3> T) , it moves as a whole in the direction of the drive. 
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Fig-EJa) shows that the magnitude of the average par- 
ticle current J is substantially larger (at the same a) 
when transport occurs through the condensed traveling 
mode, than when in the stationary mode. This can 
be understood by noticing that for well separated par- 
ticles, which are effectively not interacting, the aver- 
age drifting motion of the particles is only resisted by 
the periodic channel potential. However, when N parti- 
cles are clustered (bonded) together, then the total pin- 
ning force exerted by the channel walls on the cluster is 
/ = —^2™ =1 dU{xi)/dxi. As we show in detail below, 
the value of this net force is very sensitive to the cluster 
size, and when the length of the cluster hN is an inte- 
ger multiple of the period of the channel potential L, the 
total pinning force on the cluster vanishes, leading to a 
maximal drift velocity equal to A [2l[ . 

V. LOW FREQUENCY AC DRIVE 

In this section we discuss the behavior of the sys- 
tem when driven by an unbiased AC (square-wave) drive 
A{t) = Asgn [cos (ut)] in the low frequency limit, i.e. in 
the limit u — > 0. We focus in particular on the behavior 
of the average rectification current (J). For vanishingly 
small frequencies, (J) is obtained as the arithmetic aver- 
age of the two unidirectional currents J + and J~, with 
J ± denoting the average currents induced by the DC 
drives ±A 



A. Maximization of the rectification current 

As shown in the previous section, for constant drive 
A, increasing the pair attraction strength a, leads to the 
onset of a condensed traveling transport mode associ- 
ated with the clustering of the particles traveling in the 
direction of the drive. As the condensation sets in, the 
opposite unidirectional currents increase in magni- 
tude. However, due to the asymmetry of the channel 
potential, the condensation sets in at different values of 
a, depending on the orientation of the drive. This phe- 
nomenon is illustrated in Fig.jBJa), where the two rele- 
vant HB points are marked, for the case when h = 0.2. 
Owing to the spatial asymmetry of U(x), the depinning 
of the stationary density profile when the drive is —A, 
with current J~ , occurs at a lower value of a than when 
the drive is +A, with current J + . Therefore, when a 
is gradually increased beyond the value at the HB point 
for negative drive —A, the cycle averaged rectification 
current (J) = (1/2)(J + + J~), is negative and increases 
in absolute value, as shown in Fig.|BJb). As a is further 
increased to the value at the HB point for positive drive 
+A, the magnitude of (J) reaches a local maximum as a 
function of a. Increasing a even further results in a de- 
crease in the magnitude of (J). This occurs because the 
particles are now transported as a condensed traveling 
mode in both directions. 




FIG. 6: (Color online) (a) The unidirectional currents J 
as functions of a, calculated using the DDFT, for S = 5L, 
N = 5, A = ±1, A = 3, T = 0.5, and for h = 0.2 (solid 
line) and h — 0.16 (dashed line). Symbols "Solid squares" 
mark the corresponding Hopf bifurcation of the stationary 
density distribution for h = 0.2. Note that to left of the Hopf 
bifurcation we show the current for the stationary solutions 
and to the right for the time-periodic solution, (b) The rec- 
tification current (J) versus a, computed as the arithmetic 
mean of J + and J~ in (a) (curves) and from direct Brow- 
nian dynamics simulations with h e s = 0.16 (symbols). All 
other system parameters are as in (a). Panel (c) gives the 
unidirectional currents J obtained from Brownian dynam- 
ics simulations as functions of a for T = 0.5 (solid line) and 
T — 0.2 (dashed line). The remaining parameters are S = 10, 
N = 10, A = ±1, A = 3, and h cB = 0.16. (d) The rectifi- 
cation currents (J), computed as the arithmetic mean of the 
currents J + and J~ in (c). 



A qualitatively similar behavior of (J) is found for a 
range of different values of the particle size h. However, 
on increasing a even further, so that it is well above the 
value at the HB points, the dependence of (J) on a be- 
comes very sensitive to the value of h. For instance, in 
Fig-HJb) the rectification current attains a second mini- 
mum at around a = 6 for h = 0.2, whereas for h = 0.16 
the second minimum disappears and (J) increases mono- 
tonically as a function of a. 

To confirm the validity of the (mean field) DDFT re- 
sults, we performed Brownian dynamics computer sim- 
ulations - i.e. we numerically integrated the Langevin 
equations of motion ([3]), in order to compare with our 
DDFT results. In order to make the simulations more 
convenient to implement, we replace the hard core po- 
tential, Whx, by an equivalent, more tractable soft core 
potential, w s (xij) = e(h*/xij) , where the constants e 
and /i* can be tuned to reproduce the desired effective 
hard-core length of the potential. For fixed e and ft* the 
effective hard core length h e s of the particles becomes 
a function of a, T and, in general, also of the number 
of particles N jig]. In our simulations we set e = 0.01 
and /i* = 0.2 which corresponds to an effective hard-core 
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FIG. 7: (Color online) (a) The unidirectional currents J 
as functions of the number of particles N, in a system with 
total length S = 10L. The heavy solid and dashed curves 
correspond to the currents in the limit of very strong attrac- 
tion (a — > oo) for ft = 0.1 and h = 0.2, respectively. The 
thin dashed curves represent the currents for non-attracting 
particles (a = 0) with h — 0.1. In the inset we display the 
T — limiting values for the critical amplitude A required for 
a current to flow, as a function of hN. In (b) we display the 
rectification currents (J) obtained from the currents displayed 
in (a). 



h c ft ~ 0.16, for a — 10. Our numerical data suggests 
that the dependence of h e s on T and ./V is rather weak, 
and can therefore be neglected. 



In Fig.^b) we compare the DDFT predictions for 
h = 0.16 with the corresponding simulation results for 
h e B ~ 0.16. The first minimum in the current (J) as a 
function of a is clearly confirmed by the Brownian dy- 
namics simulation results for N — 5 particles. The sim- 
ulation results displayed in Figs-Htc) and (d) also show 
that the overall structure of (J) as a function of a does 
not change much as N is increased up to 10. Moreover, 
as the temperature is decreased from T — 0.5 down to 
T = 0.2, the maximum in the magnitude of the current 
curve, | (J) |, becomes even more pronounced, with the 
magnitude of the peak rectification current increasing by 
one order of magnitude. This effect, which is well estab- 
lished in the ratchet literature [l| , underlines the key role 
of noise in activating transport (in either direction) when 
the amplitude of the drive is smaller than both the de- 
pinning thresholds, Al and An, of the ratchet potential, 
U(x). 



B. Strong attraction limit 

In order to study the properties of the system when 
the attraction between the particles dominates over the 
thermal motion of the particles and the pinning by the 
external potential, we consider the limit a — > oo. This al- 
lows us to reduce the system of equations to a single 
equation of motion for the center of mass of the parti- 
cle condensate, y — (1/N) J^^Li x i- As noted above in 
Sec. IIVB1 the total force exerted by the channel poten- 
tial on the condensate is / = — Yli=i dU(xi)/d,Xi. If we 
assume that the pair attraction is so strong that the rods 
are closely packed together in a single condensate with 
their ends touching, the total force / can be rewritten 
as / = — J2iLidU(x + (i — l)h)/dx, where x denotes 
the coordinate of the center of the first particle in the 
file. Now, if we assume that h is small compared to the 
period L of the channel potential, then the sum can be 
replaced by an integral: 



J hN 



x+hN 



dU(y) 



dy 

U(x + hN) - U(x) 
hN 



dy 



(15) 



leading to the following effective equation of motion for 
the center of mass: 

dx _ U(x + hN) - U(x) 
~dt ~ hN 

Here, £(t) has the same statistics as in Eq. ©. To 
derive Eq. (fT6) . we use the fact that the sum of N inde- 
pendent sources of Gaussian white noise with variance 1 
is also a Gaussian noise, but with variance 1/N. 

Equation (|16p corresponds to the equation of motion 
for a single Brownian particle diffusing in the effective 
external potential V e g(x) = J[U(x + hN) — U(x)] dx, 
in contact with a thermal bath with temperature T/N . 
The first observation from Eq. (TTB1) is that for large con- 
densates, diffusion becomes negligible, so that become 
sizable only if the drive amplitude, A, overcomes the pin- 
ning force induced by the effective potential V c ft(x). For 
T = 0, this critical amplitude A is plotted in the inset 
of Fig.[?Ua), as a function of the size of the condensate 
hN. Within the shaded area, the condensate is pinned 
by the effective external potential V c ft] depinning occurs 
either to the left or to the right, depending on the drive 
orientation. Note that for hN — > 0, the right and left 
critical amplitudes coincide with the single particle de- 
pinning thresholds, Al and Ar, introduced in Sec. IIVI 
Similarly to the case for pointlike particles , selecting 
an appropriate combinations of h and N, one can achieve 
the complete locking of the condensed mode in one di- 
rection, but not in the other [2l| , which yields the upper 
bound \ A\/2 for the modulus of (J). 

Finally, using Eq. (TIE)) , we compare in FigJT] the ef- 
ficiency of the low frequency transport of strongly at- 
tracting (a — ¥ oo) and non-interacting (a = 0) particles. 
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We fix the size of the particles h and change the aver- 
age density p by changing the number of particles N in 
the system. The unidirectional currents for the conden- 
sate oscillate with N and hit the respective upper (lower) 
bound, J ± = ±\A\, for hN equal to a multiple of L. In 
the absence of particle attraction, | J* | increases mono- 
tonically with N and attain the same upper bound only 
for hN = S. The corresponding rectification currents are 
shown in FigJTfb). For certain combinations of h and N, 
the magnitude of the current of the condensate is sev- 
eral orders of magnitude larger than for non-attracting 
particles. 



VI. CONCLUDING REMARKS 

In this paper we have developed a DDFT for study- 
ing the dynamics of a file of attracting colloidal particles 
confined within a channel that exerts a periodic ratchet 
potential on the colloids. We find that the attraction 
between the colloids leads to rather rich behavior in the 
DDFT model when the particles are driven, including 
transitions from stationary to time-periodic density pro- 
files as the strength of the attraction between the par- 
ticles is increased. We also find that for strong enough 
attraction, there can be coexistence of stable stationary 
density profiles with different spatial periods and time- 
periodic density profiles, each with different values for 
the particle current J. 

These dynamical transitions in our model stem from 
the fact that the approximate free energy functional ([5]) 
on which our DDFT is based, predicts that the system 
exhibits gas-liquid phase separation for sufficiently large 
values of the ratio a/ThX. This prediction comes as a 
consequence of the mean-field approximation made in 
constructing the free energy. In reality, for a system 
containing a finite number of particles, there is no true 
phase transition. Furthermore, since the system is one- 
dimensional, there is no phase transition even in the in- 
finite sized system (i.e in the thermodynamic limit when 
N, S —> oo, with average density p = N/S remaining 
constant). In ID systems such as that studied here, as 
the attraction strength a is increased, the particles in- 
creasingly tend to gather together, but no true phase 
transition can be defined. Thus, in reality, as can be 
inferred from our Brownian dynamics simulation results, 
there are no 'sharp' transitions from the pinned to the de- 
pinned (time-periodic) state, as a is increased. Thus, we 
expect that fluctuations will round the predicted transi- 
tions. Nonetheless, as the comparison with the Brownian 
dynamics simulations show, the results from our DDFT 
do capture the main features of the system - i.e. that for 
lower values of the attraction strength a, the particles 
are uniformly distributed and that at higher values of a 
the particles gather to form a cluster, and that if the sys- 
tem length S is sufficiently long, this clustering leads to 
time-periodic currents J when the system is driven. 

In our discussions above we have pointed out that sim- 



ilarities exist between the DDFT equation (|SJ) for the 
particle density employed here and thin film equations 
that are used to model the dynamics of films and drops 
of partially wetting liquids on heterogenous solid sub- 
strates with and without additional driving forces [30| . 
The similarities result from the fact that in both cases ki- 
netic equations for conserved fields are used, and that the 
respective free energy functionals contain terms which re- 
sult in similar physical effects. For instance, the role of 
the particle-particle interactions in the present work is 
taken by wettability effects in the context of droplet dy- 
namics. The parallels between the two systems have al- 
lowed us to use the knowledge gained from studying one 
system to understand aspects of the other. In particular, 
in the present work we have drawn on the understand- 
ing of depinning mechanisms developed for thin films in 
Refs. |33l. l47j . Furthermore, Ref. [47i] indicates that one 
might encounter rich nonlinear behaviour when consid- 
ering the behaviour of attractive hard-core particles in 
wider corrugated channels, i.e., without the 'restriction' 
of single file motion. Note, however, that there are clear 
limits to the similarities: The thin film models referred to 
above do not account for any effect that is equivalent to 
the freezing instability discussed above. We believe that 
studying in detail the similarities and differences between 
DDFT and thin film models is worthwhile, as it will al- 
low for much cross fertilisation of ideas and techniques 
between the two fields. 

One of the most striking features of our system is that 
the current J depends very sensitively on the size of the 
particles h and on the total number of particles in the 
system N, particularly when the particles are strongly 
attracted to one another so that they are bound together 
to form a cluster that moves as a unit through the sys- 
tem, when there is an external drive on the system. In 
fact, the direction of travel can be completely reversed 
when the system is driven by an AC potential, simply by 
changing the number of particles in the file by one - i.e. 
adding an extra particle to a file can cause it to reverse its 
direction of motion without changing the external drive. 
This means that one can use the present system to form 
a molecular shuttle that moves back and forth between 
two docking stations, loading and unloading sing le parti- 
cles from a source to a sink docking station [21|. As the 
process repeats, a steady flux is established along the 
channel. This mechanism can be highly efficient if the 
system parameters are carefully tuned. 

The present model thus provides a useful system for 
developing a deep understanding of the behavior of 
driven macromolecular and colloidal systems occurring in 
nanoscience and biology. In particular, by using DDFT, 
which is based on a fully microscopic expression for the 
Helmholtz free energy functional ([6]) , we are able to build 
into our theory a reliable description of the correlations 
between the particles, and their influence on the dynam- 
ics of the system as a whole. 
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